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Abstract 

Complete spectra of the staggered Dirac operator ID are determined in quenched four- 
dimensional SU{2) gauge fields, and also in the presence of dynamical fermions. Periodic 
as well as antiperiodic boundary conditions are used. An attempt is made to relate the 
performance of multigrid (MG) and conjugate gradient (CG) algorithms for propagators 
with the distribution of the eigenvalues of Jp . The convergence of the CG algorithm is 
determined only by the condition number k and by the lattice size. Since /c's do not vary 
significantly when quarks become dynamic, CG convergence in unquenched fields can be 
predicted from quenched simulations. On the other hand, MG convergence is not affected 
by k but depends on the spectrum in a more subtle way. 
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1 Introduction 



Big efforts have been undertaken to find efficient multigrid (MG) methods for the computation 
of propagators in background gauge fields [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. The goal is to find an 
improved method for the solution of discretized Dirac equations. Such an improvement would 
accelerate numerical simulations of theories involving dynamical fermions considerably when 
one uses the Hybrid Monte Carlo algorithm [11]. ^ 

Although ultimately one wants to simulate theories with dynamical fermions, all the works 
on MG methods mentioned above focussed only on quenched gauge fields. However, it is 
reasonable to expect that MG methods have a chance to perform better when one considers 
"real" gauge fields which are generated in the presence of dynamical fermions. On the other 
hand one will not expect any big difference for the behavior of the conjugate gradient (CG) 
algorithm. The reasons for these two statements are as follows. The inclusion of the fermionic 
determinant in the Monte Carlo process will tend to decrease the number of (approximate) zero 
modes. This is so because configurations with less low-lying modes are more probable. MG 
methods intend to take care of the low-lying modes (which are responsible for critical slowing 
down (CSD)) on coarser grids, and the task of dealing with a reduced number of low-lying 
modes should be easier. Concerning the CG algorithm, its (asymptotic) convergence properties 
are determined by the condition number. 2 ' Since condition numbers of the (negative squared) 
massless Dirac operator are not influenced dramatically by the presence of dynamical quarks, 
one does not expect a significant consequence for the convergence behavior of CG. 

Previously other works have been done which were concerned with the role of low-lying 
modes in MG methods. The Israel group made a visual study of approximate zero modes in 
the quenched Schwinger model in the framework of their multigrid algortihm [2]. The present 
author showed that there exists an "idealized" MG algorithm which is able to eliminate CSD 
(or at least to reduce CSD strongly) both for bosonic propagators [6] and for propagators of 
staggered fermions [9] in quenched four-dimensional SU(2) gauge fields. A prerequisite for the 
success of the idealized MG method was the preservation of criticality of the Dirac operator 
under coarsening. In the present study we focus on the consequences of dynamical fermions 
for the performance of a simpler MG algorithm. This simpler algorithm proved to work in the 
quenched case, but it is unable to outperform CG [8]. Nevertheless, it is worth to see how it 
performs when dynamical fermions are present. 

This paper is organized as follows. In Sec. 2 we recall the deterministic MG method. The 
spectrum of the staggered Dirac operator is determined by means of a Lanczos procedure. 
This procedure is recalled and numerical results are presented in Sec. 3. Computations of 
propagators by CG and by MG are reported in Sec. 4, and we end with some conclusions. 

^There is a recent alternative proposal by Liischer [12] where one does not have to invert Dirac operators 
directly. 

2 )For a positive hermitean matrix the condition number equals the ratio of the largest to the smallest 
eigenvalue. 
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2 Multigrid method 



We wish to solve the squared Dirac equation 

(-^ + m 2 ) X = f (1) 

by MG, where Jp is the gauge covariant staggered Dirac operator, m is a bare mass parameter, 
and / may be a pseudofermion field, for instance. We measure physical quantities in units of a, 
where a denotes twice the lattice spacing of the staggered lattice. The reason for this is that 
free staggered fermions enjoy translational invariance only by shifts of a, not a/2 [13, 14]. 

The following MG notations will be used. The fundamental lattice is denoted by A . The 
first block lattice A 1 is obtained by coarsening with a factor of L^. Thus A 1 has Lf fewer 
sites than A (in d space-time dimensions). Restriction and interpolation operators C and A, 
respectively, are given by kernels C(x,z) and A(z,x) with z £ A , x £ A 1 . Note that C(x,z) 
and A (z, x) are N c X N c matrices in a gauge theory with N c colors. Also, C and A depend on 
the gauge field, although this is not indicated explicitly. 

In a twogrid algorithm one performs a certain number, say one, of conventional relaxation 
sweeps on A , and one obtains an approximate solution % of Eq. (1). The error e = % — % 
is unknown, but the residual r = f — (—ft 2 + rn 2 )x is computable. Error and residual are 
connected by the residual equation (—ft 2 + m 2 ) e = r. This equation is solved on A 1 where it 
reads 3 ' 

[C(-P 2 + m 2 )A]e hlock = Cr . (2) 

Solving Eq. (2) is simpler than solving the original equation because there are fewer degrees 
of freedom on A 1 . In the coarse grid correction step one replaces % by % + deblock- Then one 
performs again a relaxation sweep on A , etc. It is obvious how this twogrid algorithm can be 
generalized to a more-level MG algorithm. 

We use a blocking procedure for staggered fermions which is consistent with the lattice 
symmetries of free fermions [14]. This forces us to choose Lf, = 3. Even Lf, are not allowed. In 
four dimensions, coarsening by a factor of three reduces the number of points by 81. Therefore 
only a two-grid algorithm was implemented. This is sufficient to test the power of the MG 
method. The residual equation on the coarse grid was solved exactly by the CG algorithm. 

The averaging kernel C is chosen according to a ground-state projection definition. In the 
present work C fulfills the gauge covariant eigenvalue equation(s) [17] 

(-A NiX C*)(z,x) = \ (x)C*(z,x) (3) 

together with a normalization condition CC* = 11, and a covariance condition C(x,x) oc 11 
where x denotes the center of block x. In Eq. (3), Xo(x) is the lowest eigenvalue of — Ajv>, 

3 )This is the Galerkin choice of the residual equation on the block lattice. It assumes that e is smooth and 
can be obtained by smooth interpolation (via *4) of a suitable function e^iock on A 1 . The notion of smoothness 
in gauge fields is discussed in Ref. [8], and a more general discussion for disordered systems can be found in 
Ref. [15]; see also the recent work by Baker [16]. 
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and — Ajv> is the gauge covariant fermionic "two-link lattice Laplacian" - denned through 
Jp 2 = A + a ^ - with "Neumann boundary conditions (b.c.)". F^ is the lattice definition 
of the field strength by means of plaquette terms. Neumann b.c. means that derivative terms 
in A are omitted where one site is in block x and the other one is in a neighboring block. For 
more details we refer to Ref. [8]. 

The ground-state projection method is numerically implementable in four-dimensional non- 
abelian gauge fields [17], and since the method is gauge covariant, no gauge fixing in computa- 
tions of propagators is required. Finally we cling to a variational method where the interpolation 
operator A equals the adjoint of the restriction operator C. 

3 Spectrum of — p 2 in the presence of dynamical fermi- 
ons 

As explained in the introduction, naively one expects Jp to have less approximate zero modes in 
the presence of dynamical fermions than in the quenched case. In order to study this conjecture 
we need firstly a Hybrid Monte Carlo program, and secondly a method to determine the low- 
lying spectrum of Jp . For the generation of four-dimensional SU(2) gauge fields coupled to 
dynamical staggered fermions a FORTRAN program with vectorized CRAY code was used, 
which had been written by S. Meyer and B. Pendleton. They used this program when they 
studied the chiral transition with many fermion flavors in the SU(2) Higgs model [18]. Meyer's 
and Pendleton's program was used with four flavors of staggered fermions. The spectrum of Jp 
was determined by means of a Lanczos procedure to which we turn next. 

3.1 Lanczos procedure 

The Lanczos procedure is a technique that can be used to solve large, sparse, symmetric eigen- 
problems [19]. The method has been used in lattice field theory for a long time, see e. g. [20]. 
In course of the Lanczos iteration one generates for a given matrix A a sequence of hermitean 
tridiagonal matrices by transformations with a matrix Q whose columns are called "Lanc- 
zos vectors". These transformations have the property that the extremal eigenvalues of 
are progressively better estimates of the extremal eigenvalues of A. (In our case A = — Jp 2 ■) 
For details about the method and its convergence properties in exact arithmetic we refer to the 
literature [19]. 

In exact arithmetic, the Lanczos iteration should stop after at most n steps when A is an 
n X n matrix. In practice, however, there are severe problems [19, 21] caused by rounding errors 
and loss of orthogonality among the Lanczos vectors. As a consequence there appear so-called 
"spurious" eigenvalues [21] which are not eigenvalues of A. A clever way of identifying spurious 
eigenvalues and coping with their presence was proposed by Cullum and Willoughby [21]. In 
their algorithm one compares the eigenvalues of with the eigenvalues of a matrix T2 which 
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equals with the first row and first column deleted. If a simple eigenvalue of is also an 
eigenvalue of T2, then this eigenvalue is spurious. 

A problem which remains in the Cullum-Willoughby algorithm is that the correct multi- 
plicities of the eigenvalues of A do not emerge. This is so because a symmetric tridiagonal 
matrix does not have degenerate eigenvalues [23]. Nevertheless, in practice the matrices 
will have multiple eigenvalues which correspond to simple eigenvalues of A. They arise because 
the iteration essentially restarts itself when the numerical instabilities become too large. 4 ' 

In the present exploratory study the complete spectrum of —Jp 2 was determined. This 
was done in order to be sure about the correctness (i. e. to have no numerical uncertainties) 
in the distribution of low-lying modes. Of the several computational variants of the Lanczos 
procedure the most stable one as described in [19, Algorithm 9.2.1 and remark on p. 492] was 
implemented. Eigenvalues of the tridiagonal matrices were determined by means of the 
NAGLIB routine F02AVE [22]. 

3.2 Numerical Results in four-dimensional SU(2) gauge fields 

Let us first recall what one knows a priori about the multiplicities of the eigenvalues of —Jp 2 . 
One knows that in SU(2) gauge fields every eigenvalue is internally twofold degenerate 5 ' [6]. 
Also, — Jp 2 couples only even lattice sites to even sites, and odd sites to odd sites. Therefore 
the spectrum S of —Jp 2 equals the union of the spectra iS eve n and <S dd °f —Jp 2 restricted to 
the even and odd sublattices, respectively. One knows that iS eve n = ^odd [9]. Thus, if there 
are no further degeneracies, then on a lattice A of volume |A| there must be |A|/2 different 
eigenvalues. Their sum must equal |Tr {—Jp 2 ) = d\A\. These statements are valid for periodic 
and for antiperiodic boundary conditions. 

Concrete numerical investigations were done with the following parameters. The spectrum 
of —Jp 2 was investigated on 6 4 and 12 4 lattices. Two different kinds of boundary conditions 
(b.c.) were used. One with periodic b.c. in all directions for gauge and Fermi fields (denoted 
"per." in the following), and one with periodic b.c. for the gauge field in all directions and 
with antiperiodic b.c. on the Fermi field in time direction and periodic b.c. in spatial directions 
(denoted "antiper." in the following). The coupling (3 = 4/g 2 of the Wilson action for the 
SU(2) gauge fields was varied between 1.8 and 5.0. Quark masses m in the Hybrid Monte 
Carlo runs were chosen to be m = 0.2 and m = 0.05. These values were also used in Meyer's 
and Pendleton's work [18]; they quote m = 0.1 and m = 0.025 since they measured physical 
quantities in units of a/2. 

In Cullum's and Willoughby's Lanczos procedure j = |A| iterations were performed to 
compute (This is at least twice the number of "good" eigenvalues which can be expected.) 

4 ) From the point of view of identifying correct eigenvalues, multiple eigenvalues of numerically computed 

are welcome, because they are guaranteed to be not spurious. 

5 ) l wish to thank U.-J. Wiese for pointing out to me that this degeneracy is due to a global charge conjugation 
symmetry which is special to SU(2). 
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Figure 1: Spectrum of —Tp 2 on 12 4 lattices with periodic boundary conditions. 

The entries in the initial Lanczos vector were chosen randomly in SU(2), and the vector was 
normalized in the 2-norm. Eigenvalues of were counted as being equal when they differed 
by less than e = 10" 9 . The same e was used to identify equal (simple) eigenvalues of 
and T 2 . Proceeding in this way it turned out that in nontrivial gauge fields there seem to be 
no additional degeneracies to the ones explained above. Only for f3 = 5.0 it was impossible 
to identify |A|/2 eigenvalues whose sum equals 4|A|. For f3 = 1.8, . . . ,2.8 the method works 
perfectly. On 6 4 (12 4 ) lattices we always found 648 (10368) eigenvalues whose sum came out 
as 5184 + S 6 (82944 + S 12 ), with S 6 < 8 • 10" 9 (S 12 < 1.7 • 10" 6 ) in REAL arithmetic on a CRAY 
Y-MP. Because of the randomness of the nonvanishing off-diagonal matrix elements of —Jp 2 
this is very good evidence that the spectrum was determined exactly. 



3.2.1 Complete Spectra 

Spectra on 12 4 lattices are shown in Figs. 1 and 2. We number the different eigenvalues A& by 
k = 0, 1, 2, . . ., with Ao < Ai < . . . The data shown for finite f3 are results obtained with gauge 
fields generated by the Hybrid Monte Algorithm in the presence of dynamical fermions with a 
mass of m = 0.2. However, on the scale of the whole spectrum there is very little difference 
compared to m = 0.05 or to quenched gauge fields. The CPU time on a CRAY Y-MP is some 
50 minutes for the implemented Lanczos procedure on 12 4 lattices, and only 8 sec for 6 4 lattices. 

In Figs. 1 and 2 the data for /? = 00 (pure gauge) are not outcomes of numerical computa- 
tions, but were taken from analytical results. In the free case the eigenvectors and eigenvalues 
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12M lattice, antiperiodic b.c. in t-direction 
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Figure 2: Spectrum of —Tp 2 on 12 A lattices with antiperiodic boundary conditions for the Fermi field 
in t-direction. 

of —Jp are easily determined. For the eigenvalues one finds on a (2N) d lattice with peri- 
odic b.c. 6 ) 



A fc = 4 f sin 2 ^ with {0,1,..., [f]}. 



(4) 



The result on a (2iV) d 1 • (2r) lattice with antiperiodic b.c. in the T-direction and periodic b.c. 
in the other directions is 



d-i 



A fc = 4£smM^I +4sin 2 



V(2Av + 1) 



N J \ 2r J 

with £ {0, 1, . . . , [y] } for /j, = 1, . . . , d - 1, and k T £ {0, 1, . . . , & T>max } 
where & Tjmax = r/2 — 1 if r is even, = (r — l)/2 if r is odd. (5) 

The /? = oo values (4) and (5) are shown in Figs. 1 and 2 with their true multiplicities modulo 4 
so that they indicate the curve which the numerical data should approach for large f3. When 
one runs the Lanczos program with free fields, one obtains the eigenvalues (4) and (5) very fast 
and accurately. 

The numerical problem mentioned above which was found for f3 = 5.0 is probably due to 
the fact that the eigenvalues in the large f3 region tend to group in the clusters (4) or (5), 
respectively, where it is hard to disentangle them numerically. 



6) [f ] = N l 2 if N is even > and = ( N ~ l )l 2 if N is odd - 
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When one looks at the complete spectra, figures for 6 4 lattices look practically the same on 
the overall range as Figs. 1 and 2, when one rescales the abscissa by 16. Only for the eigenvalues 
of the free — Jp 2 this is not true. These values are collected in Table 1. 





periodic 


b.c. 






antiperiodic 


b.c. 






eigenvalue 





3 


6 


9 


12 


eigenvalue 


1 


4 


7 


10 


13 


degeneracy 


8 


64 


192 


256 


128 


degeneracy 


16 


104 


240 


224 


64 



Table 1: Spectrum of the free staggered —ID 2 on 6 4 lattices. Degeneracies are given modulo 4- 



3.2.2 Low-lying Spectra 

Now we take a closer look at the low-lying spectra of —Jp 2 . Due to renormalization effects it 
is difficult to say for which triples (A,/3,m) (+ b.c.) the results can be compared physically. 
We do not intend to do such a comparison in this exploratory study. Nevertheless, one could 
possibly discover trends, and it is instructive to look at the results as one keeps two parameters 
fixed. 

Fig. 3 shows the lowest eigenvalues (modulo the fourfold degeneracy mentioned above) of 
— Jp 2 on 6 4 lattices at (3 = 1.8. We look at individual configurations, but all Monte Carlo runs 
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m = 0.05, per. 
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Figure 3: Low-lying spectrum of —Jp 2 on 6 4 lattices at f3 = 1.8. m is the quark mass in the Hybrid 
Monte Carlo program, "antiper. " and "per. " stand for the choice of boundary conditions. The 
examples shown are configs. # 1, 2, 5, 6, 9, and 10 of Table 4- 
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were independent. In the quenched case and for m = 0.2 one is in the confined chirally broken 
phase, while for m = 0.05 one is just in the deconfined chirally symmetric phase [18]. We 
note the following. In the confined phase there is little difference in the spectra with periodic 
and with antiperiodic b.c, while in the deconfined phase (or close to the transition) there is a 
difference. At this point we also want to add the following observation. The chiral condensate 
(xx) as measured by Meyer's and Pendleton's Monte Carlo program sees no difference between 
the two kinds of b.c. (within statistical errors), in both phases. In contrast, the Polyakov loop 
is sensitive to the choice of b.c. Disregarding renormalization effects, Fig. 3 confirms the naive 
expectation that the effect of dynamical quarks is to rise the low-lying spectrum of — Jp 2 \ there 
are less approximate zero modes the lighter the fermions are. 

Figs. 4-6 show examples of the low-lying spectra of the massless operator —Jp 2 on 12 4 lattices 
for the pure gauge theory (with static quarks), and in the presence of dynamical staggered 
fermions with masses m = 0.2 and m = 0.05, respectively. Here we note again that the results 
do not depend on the choice of b.c. for smaller values of /?, while they do for larger f3. In the 
limiting case of free fields the lowest eigenvalue is always zero in case of periodic b.c. no matter 
how big the lattice is. In case of antiperiodic b.c. it equals 4 sin 2 [-7r/(2r)], see Table 2; note that 
in numerical simulations on 6 4 and 12 4 lattices the lowest eigenvalue is always much smaller 
than the free value, it is almost zero. 
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Figure 4: Low-lying spectra of the quenched operator —Tp 2 on 12^ lattices, "antiper. " and "per. " 
stand for the choice of boundary conditions. Results for f3 > 2.4 are outside the range of this plot; 
only for periodic b.c. \q comes back to zero as f3 — > oo. The examples shown are configs. # 13-16 of 
Table 4. 
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eigenvalue # 

Figure 5: Low-lying spectra of —Lp 2 on 12 A lattices in the presence of dynamical staggered fermions 
of mass 0.2. "antiper." and "per." stand for the choice of boundary conditions. (For "f3 = 2.8, per." 
only the lowest eigenvalue is visable (i.e. is < 0.001,); this point almost coincides with the one for 
"f3 = 2.4, per.".) The examples shown are configs. # 23-28 of Table 4- 



|A| 


6 4 


12 4 


18 4 


24 4 


30 4 


36 4 


72 4 


Ao 


1.0 


0.2679491 


0.1206148 


0.0681483 


0.0437048 


0.0303845 


0.0076106 



Table 2: Lowest eigenvalue of the free —Lp 2 on staggered (2L) d 1 • (2r) lattices with antiperiodic b.c. 
in t -direction, and periodic b.c. in the d — 1 spatial directions. 



3.2.3 Condition numbers 

To conclude this section let us give an idea of condition numbers. Call the two masses which 
were used in the Hybrid Monte Carlo runs mi = 0.1 and m,2 = 0.05. Denote by /c^, i = 1,2, 
the condition number of (—1ft 2 + mf), i-e- R>i = (A max + mf)/(\ m i n + mf), where A max and A m ; n 
(= Ao) denote the highest and the lowest eigenvalue of — ft 2 , respectively. Results for free fields 
are in Table 3. Table 4 gives examples in particular nontrivial configurations on 6 4 and 12 4 
lattices. (Note that it makes no sense to quote «2 for Hybrid Monte Carlo runs with mi, and 
vice versa.) The configuration numbers ("config. are referred to in Sec. 4. 

Condition numbers are (much) larger in nontrivial gauge fields than in case of free fields, 
in particular for antiperiodic b.c. This is just another manifestation that the inversion of 
(—p 2 + m 2 ) in nontrivial gauge fields is much harder than the computation of free propagators. 
Finally, note that A max in Table 4 depends in general little on the choice of boundary conditions. 
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8 



2 4 6 

eigenvalue # 

Figure 6: Low-lying spectra of —Jp 2 on 12 4 lattices in the presence of dynamical staggered fermions 
of mass 0.05. "antiper." and "per." stand for the choice of boundary conditions. (Xq = 0.011 for 
"f3 = 2.4, antiper.", so that for these parameters no point is visable on the scale of the plot.) The 
examples shown are configs. # 29-34 °f Table 4- 



|A| 


b.c. 


^min 


^max 


«i 


K2 


6 4 


per. 





12 


301 


4801 


6 4 


antiper. 


1 


13 


12.5 


13.0 


12 4 


per. 





16 


401 


6401 


12 4 


antiper. 


0.26795 


15.73205 


51.2 


58.2 



Table 3: Extremal eigenvalues of the free —p 2 , and condition numbers K{. 
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config. # 


|A| 


b.c. 


m 


P 


^min 


^max 




K2 


1 


6 4 


per. 


oo 


1.8 


1.222 


10" 


-4 


21.15 


518.1 


8066 


2 


6 4 


antiper. 


oo 


1.8 


1.725 


10" 


-4 


21.14 


527.2 


7910 


3 


6 4 


per. 


oo 


2.8 


5.668 


10" 


-2 


18.28 


189.5 


309 


4 


6 4 


antiper. 


oo 


2.8 


2.029 


10" 


-1 


18.48 


76.2 


90 


5 


6 4 


per. 


0.2 


1.8 


1.119 


10" 


-4 


20.73 


517.8 


- 


6 


6 4 


antiper. 


0.2 


1.8 


4.045 


10" 


-5 


20.76 


519.5 


- 


7 


6 4 


per. 


0.2 


2.8 


4.016 


10" 


-1 


18.25 


41.4 


— 


8 


6 4 


antiper. 


0.2 


2.8 


4.304 


10" 


-1 


18.38 


39.2 


- 


9 


6 4 


per. 


0.05 


1.8 


3.395 


10" 


-4 


20.58 


- 


7248 


10 


6 4 


antiper. 


0.05 


1.8 


4.052 


10" 


-3 


20.36 


- 


3107 


11 


6 4 


per. 


0.05 


2.8 


3.887 


10" 


-1 


18.38 


- 


47 


12 


6 4 


antiper. 


0.05 


2.8 


3.990 


10" 


-1 


18.32 


- 


46 


13 


12 4 


per. 


oo 


2.0 


1.747 


10" 


-6 


20.73 


519.2 


8286 


14 


12 4 


antiper. 


oo 


2.0 


8.574 


10" 


-7 


20.70 


518.5 


8277 


15 


12 4 


per. 


oo 


2.4 


7.556 


10" 


-5 


19.30 


482.6 


7494 


16 


12 4 


antiper. 


oo 


2.4 


1.588 


10" 


-4 


19.29 


481.3 


7255 


17 


12 4 


per. 


oo 


2.6 


4.204 


10" 


-3 


18.82 


426.7 


2807 


18 


12 4 


antiper. 


oo 


2.6 


1.916 


10" 


-2 


18.82 


318.8 


869 


19 


12 4 


per. 


oo 


2.7 


4.915 


10" 


-2 


18.66 


209.8 


361 


20 


12 4 


antiper. 


oo 


2.7 


2.621 


10" 


-2 


18.66 


282.4 


650 


21 


12 4 


per. 


oo 


2.8 


4.002 


10" 


-2 


18.47 


462.5 


434 


22 


12 4 


antiper. 


oo 


2.8 


6.198 


10" 


-2 


18.48 


181.6 


287 


23 


12 4 


per. 


0.2 


1.8 


3.075 


10" 


-6 


20.68 


514.0 


- 


24 


12 4 


antiper. 


0.2 


1.8 


6.911 


10" 


-7 


20.70 


518.5 


- 


25 


12 4 


per. 


0.2 


2.4 


8.475 


10" 


-5 


18.97 


474.2 


- 


26 


12 4 


antiper. 


0.2 


2.4 


1.602 


10" 


-5 


19.13 


479.1 


- 


27 


12 4 


per. 


0.2 


2.8 


8.250 


10" 


-5 


18.31 


457.8 


— 


28 


12 4 


antiper. 


0.2 


2.8 


3.249 


10" 


-4 


18.30 


454.8 


- 


29 


12 4 


per. 


0.05 


1.8 


4.222 


10" 


-7 


21.38 


- 


8550 


30 


12 4 


antiper. 


0.05 


1.8 


9.507 


10" 


-8 


21.55 


- 


8620 


31 


12 4 


per. 


0.05 


2.4 


1.956 


10" 


-6 


19.73 


- 


7886 


32 


12 4 


antiper. 


0.05 


2.4 


1.053 


10" 


-2 


18.93 


- 


1453 


33 


12 4 


per. 


0.05 


2.8 


5.829 


10" 


-5 


18.46 


- 


7216 


34 


12 4 


antiper. 


0.05 


2.8 


8.323 


10" 


-6 


18.88 




7527 



Table 4: Examples for the extremal eigenvalues hfi—fl 2 , and for condition numbers Ki in particular 
gauge fields on |A| lattices. The value in the column "m" gives the value of the mass of the dynamical 
fermions in the Hybrid Monte Carlo run; m = oo stands for a quenched simulation. 



4 Inversion of (— p 2 + m 2 ) 



A comprehensive summary about the computation of propagators by means of various algo- 
rithms in quenched gauge fields can be found in [8]. Here we focus on the standard CG algorithm 
[19] and the multigrid method of Sec. 2. 

4.1 Results of the conjugate gradient algorithm 

One often finds the general statement that the speed of convergence of CG depends on the 
condition number [24] . In cases where the extremal eigenvalues are well separated one can find 
"superlinear convergence", i.e. convergence at a rate that increases per iteration. More precisely 
[25], the asymptotic convergence rate of CG depends exclusively on the condition number (i.e. 
only on the extremal eigenvalues), but the form of the convergence behavior is influenced by 
the entire spectrum. If the eigenvalues are not distributed uniformly between A m ; n and A max 
(i.e. if they are clustered or there are large degeneracies), then CG converges better than the 
estimate determined by the condition number. 

In case of free fields the eigenvalues of (—ft 2 + m 2 ) are clustered, Eqs. (4) and (5), and 
the results of the computation of free propagators by CG [8] may be interpreted as a kind 
of superlinear convergence. As we saw in Sec. 3, in nontrivial gauge fields the eigenvalues are 
distributed uniformly between A m ; n and A max so that "standard" convergence must be expected. 
It is already known [8] that the inversion of (—ft 2 + m 2 ) becomes harder the more disordered 
the gauge field becomes. 

A result of the present study is that the convergence behavior of CG in nontrivial gauge fields 
is practically only determined by the condition number k, of (— Jp 2 + m 2 ), and by the lattice size; 
see Fig. 7. For configurations on a lattice of given size with the same /c, CG yields sequences of 
RMS norms of residuals which practically coincide, even if the spectra are different. This comes 
as no surprise for configurations where the spectra are almost identical (e.g configs. $ 13 and 
14, 23 and 24, 29 and 30), but it is also true in cases where there are more differences in the 
spectra (e.g. configs. $ 15 and 16, 25 and 26). However, on the other hand, as mentioned in 
Sec. 3, on the overall range of the spectra there is little difference between quenched simulations 
and simulations with dynamical fermions (of mass m = 0.2, 0.05). Therefore slight fluctuations 
in the distrubution of eigenvalues on small scales do not affect the convergence of CG. Thus if 
one wants to study the convergence of the CG algorithm one can do that with "cheap" quenched 
gauge fields, one does not have to take "expensive" unquenched configurations. 

Finally let us note that identical convergence behavior (measured by RMS norms) was found 
for m = 0.2 on 6 4 and 12 4 lattices. Only for the smaller mass m = 0.05 the RMS residual is 
reduced faster on the 6 4 lattice (Fig. 7). 
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50 100 150 200 250 300 350 400 450 500 

iteration # 

Figure 7: CG convergence of the RMS residuals in dependence on the condition number k. Curves 
1~4 are results for 12 4 lattices with k = 479,519, 1453,7886, respectively. Curves 5 and 6 are results 
on 6 4 lattices with k = 3107, 7248, respectively. The curve for convergence on a 6 4 lattice with k = 518 
coincides with curve # 2. 

4.2 Results of the twogrid algorithm 

For the inversion of (—ft 2 + m 2 ) by means of an MG method we used the twogrid algorithm 
described in Sec. 2, where the relaxation scheme on the fine grid was successive overrelaxation 
(SOR) with a relaxation parameter ou, and sweeping was done in lexicographic ordering. Ac- 
cording to the conventional MG wisdom GauB-Seidel relaxation (a; = 1) is a good smoother. 
However, from previous works [8] we know that the picture changes in nontrivial gauge fields. 
The performance of our simple variational MG method can be improved at finite gauge coupling 
by choosing co > 1. 

An obvious statement is that convergence of the MG algorithm is not determined by the 
condition number k. This is clear in the limiting case of free fields, because in pure gauges crit- 
ical slowing down is completely eliminated by MG, i.e. convergence is completely independent 
of K. 

In nontrivial gauge fields convergence of MG depends on details of the spectrum. For 
instance, configs. ft 1, 2, 5, and 6 all have the same k, for m = 0.2 (Table 4). The spectra are 
practically equal for configs. ft 1 and 2, and for configs. ft 5 and 6, and so is the MG convergence 
within each of the two groups (also as a function of a;). But MG convergence in config. ft 1 is 
distinctly different from config. ft 5. Convergence was monitored for co = 1.0, 1.2, 1.4, 1.6, 1.8, 



14 




50 100 150 200 250 300 350 400 450 500 

iteration # 

Figure 8: MG convergence of RMS residuals. The numbers at the curves have the same meaning as 
in the caption of Fig. 7. The relaxation parameter oj is 1.8 for m = 0.2, and 1.9 for m = 0.05. Here 
the curve for convergence on a 6 4 lattice with k = 518 does not coincide with curve # 2, as it does in 
Fig. 7. 

and 1.9. In conflg. $1 the best value was 1.6, while it was 1.8 in conflg. $5. MG with plain 
GauB-Seidel relaxation performed identical on configs. $1 and 5. In all cases evident inferiority 
of MG was found compared to CG, a factor of about 10 in CPU time. 

On 12 4 lattices we monitored MG convergence for the same set of a;- values mentioned above, 
and for all 12 4 configurations of Table 4. The best upvalue depends on the individual gauge 
field. Roughly speaking one obtains best convergence if one chooses co = 1.8 for m = 0.2, and 
co = 1.9 for m = 0.05. But again we could not find any significant difference in the performance 
of the MG algorithm in quenched and unquenched gauge fields. The poor performance of MG 
found earlier [8] is no feature of quenched computations. 

We conclude by giving results of MG computations in Fig. 8, where convergence is shown 
for the same configurations as in Fig. 7. We stress that we show convergence in number of 
iterations. Conversion to CPU time favors CG by another factor of 4.5. 

5 Conclusions 

The complete spectrum of the staggered Dirac operator in four-dimensional SU(2) gauge fields 
can be determined very accurately by Cullum's and Willoughby's Lanczos procedure [21], pro- 
vided the Wilson coupling (3 = A/g 2 is not too large. At finite (3 the eigenvalues of — Jp 2 
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are distributed uniformly between the lowest and the highest eigenvalue. This is so both for 
quenched simulations and for simulations with dynamical fermions. On the overall scale the 
shape of the spectrum depends little on the fermion mass. As a consequence the convergence of 
the CG algorithm is only determined by the condition number k. On a lattice of given size CG 
produces iterates whose norms depend only on k. Since k, is almost not affected by the presence 
of dynamical fermions, one can predict the convergence of CG in unquenched simulations from 
quenched simulations. 

With antiperiodic boundary conditions the lowest eigenvalue of — Jp 2 is 4 sin 2 [-7r/(2r)] which 
is not so close to zero on lattices of realistic size. However, when one introduces a nontrivial 
gauge field the lowest eigenvalue is brought very close to zero. Moreover, for intermediate values 
of f3 the spectra are practically the same for periodic and for antiperiodic boundary conditions. 

On a 6 4 lattice we found that at fixed (3 the low-lying spectrum is raised when dynamical 
fermions are introduced, and that this rise is bigger the lighter the mass of the dynamical 
quarks become. (This is a general trend, also when one passes the finite-temperature phase 
transition.) Naively this can be taken as a confirmation of the expectation that the effect of 
dynamical fermions is to suppress configurations with many approximate zero modes. However, 
one has to consider renormalization effects, which we did not intend to do in this exploratory 
study. That renormalization effects play an important role can be seen already from the results 
on 12 4 lattices. 

For the performance of the variational MG method studied here we could not find any im- 
provement when quenched gauge fields are replaced by configurations with dynamical fermions. 
We could only rediscover the previously known result [8] that there will be a breakeven point 
in lattice sizes after which MG will outperform CG. This is so because in the limiting case 
f3 — y oo critical slowing down is completely eliminated. However, we cannot judge how big the 
lattices have to be in order to reach the breakeven point. It is reasonable to believe that in 
principle the presence of dynamical fermions will affect the performance of MG algorithms in 
a positive way. We think that the main reason for seeing no improvement is that the notion of 
"Laplace smoothness" [14] which stands behind the definition (3) is inappropriate for staggered 
fermions. One should rather use the "Dirac notion of smoothness" [14, 10] which is in the spirit 
of the discussion in Ref. [15, 16]. Possibly with Baker's algorithm [16] one can observe that the 
presence of dynamical fermions simplifies the task for MG algorithms. 
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